{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "nbsphinx": "hidden"
   },
   "source": [
    "This notebook is part of the `kikuchipy` documentation https://kikuchipy.org.\n",
    "Links to the documentation won't work from the notebook."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Feature maps\n",
    "\n",
    "This section details methods for extracting information from pattern\n",
    "intensities, called feature maps (for lack of a better description).\n",
    "\n",
    "Let's import the necessary libraries and a Nickel EBSD test data set:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Exchange inline for notebook or qt5 (from pyqt) for interactive plotting\n",
    "%matplotlib inline\n",
    "\n",
    "import hyperspy.api as hs\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import kikuchipy as kp\n",
    "\n",
    "\n",
    "# Use kp.load(\"data.h5\") to load your own data\n",
    "s = kp.data.nickel_ebsd_large(allow_download=True)  # External download\n",
    "s"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Image quality\n",
    "\n",
    "The image quality metric $Q$ presented by Krieger Lassen\n",
    "<cite data-cite=\"lassen1994automated\"> Lassen (1994)</cite>\n",
    "can be calculated for an [EBSD](reference.rst#kikuchipy.signals.EBSD) signal\n",
    "with\n",
    "[get_image_quality()](reference.rst#kikuchipy.signals.EBSD.get_image_quality),\n",
    "or, for a single pattern (`numpy.ndarray`), with\n",
    "[get_image_quality()](reference.rst#kikuchipy.pattern.get_image_quality).\n",
    "Following the notation in\n",
    "<cite data-cite=\"marquardt2017quantitative\">Marquardt et al. (2017)</cite>, it\n",
    "is given by\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "Q &= 1 - \\frac{J}{J_{\\mathrm{res}}w_{\\mathrm{tot}}},\\\\\n",
    "J &= \\sum_{h = -N/2}^{N/2} \\sum_{k = -N/2}^{N/2} w(h, k)\n",
    "\\left|\\mathbf{q}\\right|^2,\\\\\n",
    "J_{\\mathrm{res}} &= \\frac{1}{N^2} \\sum_{h = -N/2}^{N/2}\n",
    "\\sum_{k = -N/2}^{N/2} \\left|\\mathbf{q}\\right|^2,\\\\\n",
    "w_{\\mathrm{tot}} &= \\sum_{h = -N/2}^{N/2} \\sum_{k = -N/2}^{N/2} w(h, k).\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "The function $w(h, k)$ is the Fast Fourier Transform (FFT) power spectrum of the\n",
    "EBSD pattern, and the vectors $\\mathbf{q}$ are the frequency vectors with\n",
    "components $(h, k)$. The sharper the Kikuchi bands, the greater the high\n",
    "frequency content of the power spectrum, and thus the closer $Q$ will be to\n",
    "unity.\n",
    "\n",
    "Since we want to use the image quality metric to determine how pronounced the\n",
    "Kikuchi bands in our patterns are, we first remove the static and dynamic\n",
    "background:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "s.remove_static_background()\n",
    "s.remove_dynamic_background()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To visualize parts of the computation, we compute the power spectrum of a\n",
    "pattern in the Nickel EBSD data set and the frequency vectors, shift the\n",
    "zero-frequency components to the centre, and plot them:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p = s.inav[20, 11].data\n",
    "p_fft = kp.pattern.fft(p, shift=True)\n",
    "q = kp.pattern.fft_frequency_vectors(shape=p.shape)\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(22, 6), ncols=3)\n",
    "title_kwargs = dict(fontsize=22, pad=15)\n",
    "fig.subplots_adjust(wspace=0.05)\n",
    "\n",
    "im0 = ax[0].imshow(p, cmap=\"gray\")\n",
    "ax[0].set_title(\"Pattern\", **title_kwargs)\n",
    "fig.colorbar(im0, ax=ax[0])\n",
    "\n",
    "im1 = ax[1].imshow(np.log(kp.pattern.fft_spectrum(p_fft)), cmap=\"gray\")\n",
    "ax[1].set_title(\"Log of shifted power spectrum of FFT\", **title_kwargs)\n",
    "fig.colorbar(im1, ax=ax[1])\n",
    "\n",
    "im2 = ax[2].imshow(np.fft.fftshift(q), cmap=\"gray\")\n",
    "ax[2].set_title(r\"Shifted frequency vectors $q$\", **title_kwargs)\n",
    "_ = fig.colorbar(im2, ax=ax[2])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we don't want the EBSD patterns to be\n",
    "[zero-mean normalized](pattern_processing.ipynb#normalize-intensity) before\n",
    "computing $Q$, we must pass `get_image_quality(normalize=False)`.\n",
    "\n",
    "Let's compute the image quality $Q$ and plot it for the entire data set:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "iq = s.get_image_quality(normalize=True)  # Default"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(10, 6))\n",
    "plt.imshow(iq, cmap=\"gray\")\n",
    "plt.colorbar(label=r\"Image quality, $Q$\", pad=0.01)\n",
    "_ = plt.axis(\"off\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we want to use this map to navigate around in when plotting patterns, we can\n",
    "easily do that as explained in the\n",
    "[visualizing patterns](visualizing_patterns.rst) guide."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Average dot product\n",
    "\n",
    "The average dot product, or normalized cross-correlation when centering each\n",
    "pattern's intensity about zero and normalizing the intensities to a standard\n",
    "deviation $\\sigma$ of 1 (which is the default behaviour), between each pattern\n",
    "and their four nearest neighbours, can be obtained for an\n",
    "[EBSD](reference.rst#kikuchipy.signals.EBSD) signal with\n",
    "[get_average_neighbour_dot_product_map()](reference.rst#kikuchipy.signals.EBSD.get_average_neighbour_dot_product_map)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "adp = s.get_average_neighbour_dot_product_map()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbsphinx-thumbnail": {
     "tooltip": "Extraction of information from pattern intensities"
    },
    "tags": [
     "nbsphinx-thumbnail"
    ]
   },
   "outputs": [],
   "source": [
    "plt.figure(figsize=(10, 6))\n",
    "plt.imshow(adp, cmap=\"gray\")\n",
    "plt.colorbar(label=\"Average dot product\", pad=0.01)\n",
    "_ = plt.axis(\"off\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The map displays how similar each pattern is to its neighbours. Grain\n",
    "boundaries, and some scratches on the sample, can be clearly seen as pixels with\n",
    "a lower value, signifying that they are more dissimilar to their neighbouring\n",
    "pixels, as the ones within grains where the neighbour pixel similarity is high.\n",
    "\n",
    "The map above was created by averaging the dot product matrix per map point,\n",
    "created by calculating the dot product between each pattern and their four\n",
    "nearest neighbours, which can be seen in the black spots (uneven sample surface)\n",
    "in the left grains"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "w1 = kp.filters.Window()\n",
    "w1.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We could instead average with e.g. the eight nearest neighbours"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "w2 = kp.filters.Window(window=\"rectangular\", shape=(3, 3))\n",
    "w2.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "adp2 = s.get_average_neighbour_dot_product_map(window=w2)\n",
    "\n",
    "plt.figure(figsize=(10, 6))\n",
    "plt.imshow(adp2, cmap=\"gray\")\n",
    "plt.colorbar(label=\"Average dot product\", pad=0.01)\n",
    "_ = plt.axis(\"off\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the window coefficients must be integers.\n",
    "\n",
    "We can also control whether pattern intensities should be centered about zero\n",
    "and/or whether they should be normalized prior to calculating the dot products\n",
    "by passing `zero_mean=False` and/or `normalize=False`. These are `True` by\n",
    "default. The data type of the output map, 32-bit floating point by default,\n",
    "can be set by passing e.g. `dtype_out=np.float64`.\n",
    "\n",
    "We can obtain the dot product matrices per map point, that is the matrices\n",
    "before they are averaged, with\n",
    "[get_neighbour_dot_product_matrices()](reference.rst#kikuchipy.signals.EBSD.get_neighbour_dot_product_matrices).\n",
    "Let's see similar a pattern on a grain boundary in map point (x, y) = (50, 19)\n",
    "is to all its nearest neighbour in a (5, 5) window centered on that point"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "w3 = kp.filters.Window(\"rectangular\", shape=(5, 5))\n",
    "dp_matrices = s.get_neighbour_dot_product_matrices(window=w3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x, y = (50, 19)\n",
    "s_dp_matrices = hs.signals.Signal2D(dp_matrices)\n",
    "s_dp_matrices.inav[x, y].plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see that the pattern is more similar to the patterns up to the right,\n",
    "while it is quite dissimilar to the patterns to the lower left. Let's visualize\n",
    "this more clearly, as is done e.g. in Fig. 1 by\n",
    "<cite data-cite=\"brewick2019nlpar\">Brewick et al. (2019)</cite>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_n, x_n = w3.n_neighbours\n",
    "\n",
    "s2 = s.inav[x - x_n:x + x_n + 1, y - y_n:y + y_n + 1].deepcopy()\n",
    "s2.rescale_intensity(percentiles=(0.5, 99.5))  # Stretch the contrast a bit\n",
    "s3 = s2 * s_dp_matrices.inav[x, y].T  # Signals must have same navigation shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "_ = hs.plot.plot_images(\n",
    "    images=s3,\n",
    "    per_row=5,\n",
    "    label=None,\n",
    "    suptitle=None,\n",
    "    axes_decor=None,\n",
    "    colorbar=None,\n",
    "    vmin=int(s3.data.min()),\n",
    "    vmax=int(s3.data.max()),\n",
    "    padding=dict(wspace=0, hspace=-0.05),\n",
    "    fig=plt.figure(figsize=(10, 10))\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we can pass this dot product matrix directly to\n",
    "[get_average_neighbour_dot_product_map()](reference.rst#kikuchipy.signals.EBSD.get_average_neighbour_dot_product_map)\n",
    "via the `dp_matrices` parameter to obtain the average dot product map from these\n",
    "matrices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "adp3 = s.get_average_neighbour_dot_product_map(dp_matrices=dp_matrices)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's plot this and highlight the location of the pattern on the grain boundary\n",
    "above with a red circle"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(10, 6))\n",
    "plt.imshow(adp3, cmap=\"gray\")\n",
    "plt.colorbar(label=\"Average dot product\", pad=0.01)\n",
    "plt.scatter(x=x, y=y, marker=\"o\", c=\"r\", s=50)\n",
    "_ = plt.axis(\"off\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
